/*
 * $Id: spm_resels_vol.c 938 2007-10-12 19:09:31Z john $
 * John Ashburner
 */
 
/*
See Worsley et al. (1996), Human Brain Mapping 4:58-73 for a description
of what it does.
*/

#include <math.h>
#include "mex.h"
#include "spm_mapping.h"

static void resel_fun(int *curr, int *prev, /* current and previous planes */
	int m, int n, /* image dimensions */
	int *P,   /* # points */
	int E[3], /* # edges  */
	int F[4], /* # faces  */
	int *C)   /* # cubes  */
{
	int p=0,ex=0,ey=0,ez=0,fxy=0,fxz=0,fyz=0,c=0;
	int i, j;
	for(i=1; i<m; i++)
		for(j=1; j<n; j++)
		{
			int o = i+j*m;

			if (curr[o])
			{
				p ++;

				/* A simple way of doing it... 
				if (curr[o-1]) ex ++;
				if (curr[o-m]) ey ++;
				if (prev[o  ]) ez ++;
				if (curr[o-1] && curr[o-m] && curr[o-m-1]) fxy ++;
				if (curr[o-1] && prev[o  ] && prev[o  -1]) fxz ++;
				if (curr[o-m] && prev[o  ] && prev[o  -m]) fyz ++;
				if (curr[o-1] && curr[o-m] && curr[o-m-1] && prev[o]
				 && prev[o-1] && prev[o-m] && prev[o-m-1]) c ++;
				*/

				/* It could be optimized much further with "if else" statements
				   but it was beginning to hurt my head */
				if (curr[o-1])
				{
					ex ++;
					if (curr[o-m] && curr[o-m-1])
					{
						fxy ++;
						if (prev[o] && prev[o-1] && prev[o-m] && prev[o-m-1]) c ++;
					}
				}
				if (curr[o-m])
				{
					ey ++;
					if (prev[o  ] && prev[o  -m]) fyz ++;
				}
				if (prev[o  ])
				{
					ez ++;
					if (curr[o-1] && prev[o  -1]) fxz ++;
				}
			}
		}
	*P   += p;
	E[0] += ex;
	E[1] += ey;
	E[2] += ez;
	F[0] += fxy;
	F[1] += fxz;
	F[2] += fyz;
	*C   += c;
}


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	MAPTYPE *map, *get_maps();
	int m,n,k,i;
	int E[3], F[3], P, C;
	double *R, r[3], *img;
	int *curr, *prev, *tmpp;

	if (nrhs != 2 || nlhs > 1)
	{
		mexErrMsgTxt("Incorrect usage.");
	}

	map = get_maps(prhs[0], &n);
	if (n!=1)
	{
		free_maps(map, n);
		mexErrMsgTxt("Bad image handle dimensions.");
	}

	if (!mxIsNumeric(prhs[1]) || mxIsComplex(prhs[1]) ||
	    !mxIsDouble(prhs[1]))
	{
		free_maps(map, 1);
		mexErrMsgTxt("Second argument must be numeric, real, full and double.");

	}

	if (mxGetM(prhs[1])*mxGetN(prhs[1]) != 3)
	{
		free_maps(map, n);
		mexErrMsgTxt("Second argument must contain three elements.");
	}

	r[0] = 1.0/mxGetPr(prhs[1])[0];
	r[1] = 1.0/mxGetPr(prhs[1])[1];
	r[2] = 1.0/mxGetPr(prhs[1])[2];

	plhs[0] = mxCreateDoubleMatrix(4,1,mxREAL);
	R       = mxGetPr(plhs[0]);

	m = map->dim[0];
	n = map->dim[1];
	k = map->dim[2];

	curr = (int *)mxCalloc((m+1)*(n+1),sizeof(int)); /* current plane */
	prev = (int *)mxCalloc((m+1)*(n+1),sizeof(int)); /* previous plane */
	img  = (double *)mxCalloc((m+1)*(n+1),sizeof(double));

	P = C = E[0] = E[1] = E[2] = F[0] = F[1] = F[2] = 0;
	for(i=0; i<k; i++)
	{
		int i1, j1;
		static double mat[16] = {
			1.0, 0.0, 0.0, 0.0,
			0.0, 1.0, 0.0, 0.0,
			0.0, 0.0, 1.0, 0.0,
			0.0, 0.0, 0.0, 1.0};
		mat[14] = i+1.0;
		slice(mat, img, m, n, map, 0, 0);
		for(i1=0;i1<m; i1++)
			for(j1=0; j1<n; j1++)
				if (mxIsFinite(img[i1+m*j1]) && img[i1+m*j1])
					curr[i1+1+(m+1)*(j1+1)] = 1;
				else
					curr[i1+1+(m+1)*(j1+1)] = 0;

		/* count edges, faces etc */
		resel_fun(curr,prev,m+1,n+1, &P, E, F, &C);

		/* make current plane previous */
		tmpp = prev; prev = curr; curr = tmpp;
	}

	(void)mxFree((char *)curr);
	(void)mxFree((char *)prev);
	(void)mxFree((char *)img);
	free_maps(map, 1);

	R[0] = P - (E[0]+E[1]+E[2])+(F[0]+F[1]+F[2])-C;
	R[1] = (E[0]-F[0]-F[1]+C)*r[0] + (E[1]-F[0]-F[2]+C)*r[1] + (E[2]-F[1]-F[2]+C)*r[2];
	R[2] = (F[0]-C)*r[0]*r[1] + (F[1]-C)*r[0]*r[2] + (F[2]-C)*r[1]*r[2];
	R[3] = C*r[0]*r[1]*r[2];
}
